home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / oper_sys / weyl / weyl_lsp.lha / gfp.lisp < prev    next >
Encoding:
Text File  |  1991-10-24  |  17.0 KB  |  505 lines
  1. ;;; -*- Mode:Lisp; Package:Weyli; Syntax:Common-Lisp; Base:10; Lowercase:T -*-
  2. ;;; ===========================================================================
  3. ;;;                  GF(p)
  4. ;;; ===========================================================================
  5. ;;; (c) Copyright 1989, 1991 Cornell University
  6.  
  7. ;;; $Id: gfp.lisp,v 2.11 1991/10/24 19:19:36 rz Exp $
  8.  
  9. (in-package "WEYLI")
  10.  
  11. (define-domain-element-classes GFp GFp-element)
  12.  
  13. (defmethod number-of-elements ((domain GFp))
  14.   (characteristic domain))
  15.  
  16. (defmethod number-of-elements ((domain GFq))
  17.   (expt (characteristic domain) (field-degree domain)))
  18.  
  19. (defmethod make-GFp-domain ((characteristic integer) (degree integer))
  20.   (cond ((= degree 1)
  21.      (make-instance 'gfp :characteristic characteristic))
  22.      (t (error "Can't do GF(~D^~D) yet" characteristic degree)
  23.        ;; This is where GFq domains are to be defined.
  24.        )))
  25.  
  26. (defmethod print-object ((d GFp) stream)
  27.   #+Genera
  28.   (format stream "~'bGF~(~D)" (characteristic d))
  29.   #-Genera
  30.   (format stream "GF(~D)" (characteristic d)))
  31.  
  32. (defmethod make-element ((domain GFp) (value integer) &rest ignore)
  33.   (declare (ignore ignore))
  34.   (let ((modulus (characteristic domain)))
  35.     (make-instance 'GFp-element
  36.            :domain domain
  37.            :value (reduce-modulo-integer value modulus))))
  38.  
  39. ;; Could have more error checking
  40. (defmethod weyl::make-element ((domain GFp) (value integer) &rest ignore)
  41.   (declare (ignore ignore))
  42.   (make-element domain value))
  43.  
  44. (defmethod print-object ((x GFp-element) stream)
  45.   (with-slots (value domain) x
  46.     (format stream "~D(~D)" value (characteristic domain))))
  47.  
  48. (defmethod = ((x GFp-element) (y GFp-element))
  49.   (with-slots ((v1 value) (d1 domain)) x
  50.     (with-slots ((v2 value) (d2 domain)) y
  51.       (and (eq d1 d2) (eql v1 v2)))))
  52.  
  53. (defmethod 0? ((x GFp-element))
  54.   (with-slots (value) x
  55.     (lisp:zerop value)))
  56.  
  57. (defmethod 1? ((x GFp-element))
  58.   (with-slots (value) x
  59.     (eql value 1)))
  60.  
  61. ;; The following three methods make finite fields behave like quotient fields
  62.  
  63. (defmethod make-quotient-element ((domain GFp) (a GFp-element) (b GFp-element))
  64.   (unless (eql domain (domain-of a))
  65.     (error "~S should be an element of ~S" a domain))
  66.   (unless (eql domain (domain-of b))
  67.     (error "~S should be an element of ~S" b domain))
  68.   (with-slots ((v1 value)) a
  69.     (with-slots ((v2 value)) b
  70.       (with-slots (characteristic) domain
  71.     (make-element domain
  72.       (lisp:* v1 (compute-inverse v2 characteristic)))))))
  73.  
  74. (defmethod numerator ((a GFp-element))
  75.   a)
  76.  
  77. (defmethod denominator ((a GFp-element))
  78.   (make-element (domain-of a) 1))
  79.  
  80. (defmethod minus ((x GFp-element))
  81.   (with-slots (value domain) x
  82.     (with-slots (characteristic) domain
  83.       (if (eql 2 characteristic) x
  84.       (make-element domain (lisp:- characteristic value))))))
  85.  
  86. ;;; There is no such thing as a negative number in finite fields.
  87. (defmethod minus? ((x GFp-element))
  88.   nil)
  89.  
  90. (defmethod plus? ((x GFp-element))
  91.   (not (0? x)))
  92.  
  93. (defmethod-binary plus GFp-element  (a b)
  94.   (make-element domain (lisp:+ (gfp-value a) (gfp-value b))))
  95.  
  96. (defmethod difference ((a GFp-element) (b GFp-element))
  97.   (with-slots ((v1 value) (d1 domain)) a
  98.     (with-slots ((v2 value) (d2 domain)) b
  99.       (cond ((not (eql d1 d2))
  100.          (error "~S and ~S are not from the same domain" a b))
  101.         (t (make-element d1 (lisp:- v1 v2)))))))
  102.  
  103. (defmethod times ((a GFp-element) (b GFp-element))
  104.   (with-slots ((v1 value) (d1 domain)) a
  105.     (with-slots ((v2 value) (d2 domain)) b
  106.       (cond ((not (eql d1 d2))
  107.          (error "~S and ~S are not from the same domain" a b))
  108.         (t (make-element d1 (lisp:* v1 v2)))))))
  109.  
  110. (defmethod plus ((a GFp-element) (b integer))
  111.   (with-slots ((v1 value) (d1 domain)) a
  112.     (make-element d1
  113.       (lisp:+ v1 (reduce-modulo-integer b (characteristic d1))))))
  114.  
  115. (defmethod plus ((a integer) (b GFp-element))
  116.   (with-slots ((v1 value) (d1 domain)) b
  117.     (make-element d1
  118.       (lisp:+ (reduce-modulo-integer a (characteristic d1)) v1))))
  119.  
  120. (defmethod difference ((a GFp-element) (b integer))
  121.   (with-slots ((v1 value) (d1 domain)) a
  122.     (make-element d1
  123.       (lisp:- v1 (reduce-modulo-integer b (characteristic d1))))))
  124.  
  125. (defmethod difference ((a integer) (b GFp-element))
  126.   (with-slots ((v1 value) (d1 domain)) b
  127.     (make-element d1
  128.       (lisp:- (reduce-modulo-integer a (characteristic d1)) v1))))
  129.  
  130. (defmethod times ((a GFp-element) (b integer))
  131.   (with-slots ((v1 value) (d1 domain)) a
  132.     (make-element d1
  133.       (lisp:* v1 (reduce-modulo-integer b (characteristic d1))))))
  134.  
  135. (defmethod times ((a integer) (b GFp-element))
  136.   (with-slots ((v1 value) (d1 domain)) b
  137.     (make-element d1
  138.       (lisp:* (reduce-modulo-integer a (characteristic d1)) v1))))
  139.  
  140. ;; Takes the inverse of an integer N mod P.  Solve N*X + P*Y = 1.  N
  141. ;; is guaranteed to be less than P, since in the case where P is a
  142. ;; fixnum, N is also assumed to be one.
  143.  
  144. (defmethod recip ((x GFp-element))
  145.   (with-slots (value domain) x
  146.     (with-slots (characteristic) domain
  147.       (make-element domain (reduce-modulo-integer
  148.                 (compute-inverse value characteristic)
  149.                 characteristic)))))
  150.  
  151. (defun compute-inverse (value modulus)
  152.   (let ((a1 modulus)
  153.     (a2 (if (lisp:< value 0) (lisp:+ value modulus) value))
  154.     (y1 0)
  155.     (y2 1)
  156.     q)
  157.     (loop
  158.       (if (eql a2 1) (return (values y2 y1)))
  159.       (if (lisp:zerop a2)
  160.       (error "Inverse of zero divisor -- ~d modulo ~d"
  161.          value modulus))
  162.       (setq q (truncate a1 a2))
  163.       (psetq a1 a2 a2 (lisp:- a1 (lisp:* a2 q)))
  164.       (psetq y1 y2 y2 (lisp:- y1 (lisp:* y2 q))))))
  165.     
  166. (defmethod expt ((x GFp-element) (e integer))
  167.   (with-slots (value domain) x
  168.     (cond ((eql 1 value) x)
  169.       ((lisp:minusp e)
  170.        (error "Raising ~D to a negative power ~D" x e))
  171.       (t (make-element domain
  172.            (expt-modulo-integer value e (characteristic domain)))))))
  173.  
  174. (defmethod quotient ((a GFp-element) (b GFp-element)) 
  175.   (with-slots ((v1 value) (d1 domain)) a
  176.     (with-slots ((v2 value) (d2 domain)) b
  177.       (cond ((eq d1 d2)
  178.          (with-slots (characteristic) d1
  179.            (make-element d1
  180.          (lisp:* v1 (compute-inverse v2 characteristic)))))
  181.         (t (error "Taking the quotient of elements of ~
  182.                different fields: ~S, ~S"
  183.               a b))))))
  184.  
  185. (defmethod remainder ((a GFp-element) (b GFp-element))
  186.   (error "Computing the remainder of ~D by ~D"
  187.      a b))
  188.  
  189. (defmethod gcd ((a GFp-element) (b GFp-element))
  190.   (with-slots ((v1 value) (d1 domain)) a
  191.     (with-slots ((v2 value) (d2 domain)) b
  192.       (cond ((eq d1 d2) (make-element d1 1))
  193.         (t (error "Taking the GCD of elements of different fields: ~S, ~S"
  194.               a b))))))
  195.  
  196. (defmethod lcm ((a GFp-element) (b GFp-element))
  197.   (with-slots ((v1 value) (d1 domain)) a
  198.     (with-slots ((v2 value) (d2 domain)) b
  199.       (cond ((eq d1 d2) (make-element d1 1))
  200.         (t (error "Taking the LCM of elements of different fields: ~S, ~S"
  201.               a b))))))
  202.  
  203. (defmethod random ((domain GFp))
  204.   (make-element domain (lisp:random (characteristic domain))))
  205.  
  206. (defmethod multiplicative-order ((a GFp-element))
  207.   (with-slots (value domain) a
  208.     (with-slots ((p  characteristic)) domain
  209.       (cond ((not (eql 1 (lisp:gcd value p)))
  210.          *positive-infinity*)
  211.         ((let ((group-order (totient p)))
  212.            (do ((factors (factor group-order)
  213.                  (rest factors))
  214.             (order group-order))
  215.            ((null factors)
  216.             order)
  217.          (do ((i 0 (lisp:1+ i)))
  218.              ((lisp:= i (cdar factors)))
  219.            (setq order (lisp:/ order (caar factors)))
  220.            (when (not (eql 1 (expt-modulo-integer value order p)))
  221.              (setq order (lisp:* order (caar factors)))
  222.              (return t))))))))))
  223.  
  224. ;; GF(2^n)
  225. (defvar *GF2-irreducible-polynomials*
  226.     '(#O7 #O13 #O23 #O45 #O103 #O211 #O435 #O1021 #O2011 #O4005 #O10123
  227.       #O20033 #O42103 #O100003 #O210013))
  228.  
  229. (defmethod make-GFp-domain ((characteristic (eql 2)) (degree integer))
  230.   (cond ((= degree 1)
  231.      (make-instance 'gfp :characteristic characteristic))
  232.     ((< degree (+ (length *GF2-irreducible-polynomials*) 2))
  233.      (let* ((mask (ash 1 degree))
  234.         (field (1- mask))
  235.         (min-poly (logand (nth (- degree 2) *GF2-irreducible-polynomials*)
  236.                   field)))
  237.        (make-instance 'GF2^n 
  238.              :degree degree
  239.          :reduction-table
  240.           (loop for i below degree
  241.             for x^n = min-poly then (ash x^n 1)
  242.             collect
  243.              (if (lisp:zerop (logand mask x^n)) x^n
  244.              (setq x^n (logxor (logand field x^n) min-poly))))
  245.           :characteristic characteristic)))
  246.     (t (error "Table doesn't go far enough: 2^~D" degree))))
  247.  
  248. (defmethod print-object ((domain GF2^n) stream)
  249.   #+Genera
  250.   (format stream "~'bGF~(2^~D)" (field-degree domain))
  251.   #-Genera
  252.   (format stream "GF(2^~D)" (field-degree domain)))
  253.  
  254. (defclass GF2^n-element (GFp-element) 
  255.      ())
  256.  
  257. (defmethod print-object ((elt GF2^n-element) stream)
  258.   (format stream "~V,'0B(2^~D)"
  259.       (field-degree (domain-of elt)) (GFp-value elt)
  260.     (field-degree (domain-of elt))))
  261.  
  262. (defmethod make-element ((domain GF2^n) (value integer) &rest ignore)
  263.   (declare (ignore ignore))
  264.   (make-instance 'GF2^N-element
  265.          :domain domain
  266.          :value (logand (1- (ash 1 (field-degree domain))) value)))
  267.  
  268. ;; Could have more error checking
  269. (defmethod weyl::make-element ((domain GF2^n) (value integer) &rest ignore)
  270.   (declare (ignore ignore))
  271.   (make-element domain value))
  272.  
  273. (defmethod multiplicative-order ((a GF2^n-element))
  274.   (let ((group-size (1- (number-of-elements (domain-of a)))))
  275.     (loop for order in (all-divisors group-size)
  276.       do (when (1? (expt a order))
  277.            (return order)))))
  278.  
  279. (defmethod-binary plus GF2^n-element (a b)
  280.   (make-element domain (logxor (gfp-value a) (gfp-value b))))
  281.  
  282. (defmethod-binary times GF2^n-element (a b)
  283.   (let ((x (Gfp-value a))
  284.     (y (GFp-value b))
  285.     (degree (field-degree domain))
  286.     (acc 0) answer)
  287.     (loop while (not (lisp:zerop y)) do
  288.       (when (not (lisp:zerop (lisp:logand 1 y)))
  289.     (setq acc (lisp:logxor acc x)))
  290.       (setq x (lisp:ash x 1))
  291.       (setq y (lisp:ash y -1)))
  292.     (setq answer (lisp:logand (lisp:1- (lisp:ash 1 degree)) acc))
  293.     (loop for hi-bits = (lisp:ash acc (lisp:- degree))
  294.         then (lisp:ash hi-bits -1)
  295.       for poly in (GFp-reduction-table domain)
  296.       while (not (lisp:zerop hi-bits))
  297.       do (unless (lisp:zerop (lisp:logand 1 hi-bits))
  298.            (setq answer (lisp:logxor answer poly))))
  299.     (make-instance 'GF2^N-element :domain domain :value answer)))    
  300.  
  301. (defmethod expt ((base GF2^n-element) (expt integer))
  302.   (%funcall (repeated-squaring #'times (make-element (domain-of base) 1))
  303.         base expt))
  304.  
  305. (defmethod recip ((x GF2^n-element))
  306.   (let ((domain (domain-of x)))
  307.     (expt x (lisp:- (lisp:expt 2 (field-degree domain)) 2))))
  308.  
  309. (defmethod-binary quotient GF2^n-element (x y)
  310.   (* x (recip y)))
  311.  
  312. ;; GF(m)
  313.  
  314. ;; This domain is the union of all Z/mZ for all m.
  315.  
  316. (define-domain-element-classes GFm GFm-element)
  317.  
  318. (defmethod make-element ((domain GFm) value &rest rest)
  319.   (let ((modulus (first rest)))
  320.     (make-instance 'GFm-element
  321.            :domain domain
  322.            :value (reduce-modulo-integer value modulus)
  323.            :modulus modulus)))
  324.  
  325. ;; Could have more error checking
  326. (defmethod weyl::make-element ((domain GFm) value &rest rest)
  327.   (%apply #'make-element domain value rest))
  328.  
  329. (defmethod print-object ((x GFm-element) stream)
  330.   (with-slots (value modulus) x
  331.     (format stream "~D(~D)" value modulus)))
  332.  
  333. (defmethod = ((x GFm-element) (y GFm-element))
  334.   (with-slots ((v1 value) (m1 modulus) (d1 domain)) x
  335.     (with-slots ((v2 value) (m2 modulus) (d2 domain)) y
  336.       (and (eq d1 d2) (eql v1 v2) (eql m1 m2)))))
  337.  
  338. (defmethod 0? ((x GFm-element))
  339.   (with-slots (value) x
  340.     (lisp:zerop value)))
  341.  
  342. (defmethod 1? ((x GFm-element))
  343.   (with-slots (value) x
  344.     (eql value 1)))
  345.  
  346. (defmethod minus ((x GFm-element))
  347.   (with-slots (value modulus domain) x
  348.     (if (eql 2 modulus) x
  349.     (make-element domain (lisp:- modulus value) modulus))))
  350.  
  351. ;;; There is no such thing as a negative number in finite fields.
  352. (defmethod minus? ((x GFm-element))
  353.   nil)
  354.  
  355. (defmethod plus? ((x GFm-element))
  356.   (not (0? x)))
  357.  
  358. (defmethod plus ((a GFm-element) (b GFm-element))
  359.   (with-slots ((v1 value) (m1 modulus) (d1 domain)) a
  360.     (with-slots ((v2 value) (m2 modulus) (d2 domain)) b
  361.       (cond ((not (eql d1 d2))
  362.          (error "~S and ~S are not from the same domain" a b))
  363.         ((eql m1 m2)
  364.          (make-element d1 (lisp:+ v1 v2) m1))
  365.         (t (make-element d1 (lisp:+ v1 v2) (lisp:gcd m1 m2)))))))
  366.  
  367. (defmethod difference ((a GFm-element) (b GFm-element))
  368.   (with-slots ((v1 value) (m1 modulus) (d1 domain)) a
  369.     (with-slots ((v2 value) (m2 modulus) (d2 domain)) b
  370.       (cond ((not (eql d1 d2))
  371.          (error "~S and ~S are not from the same domain" a b))
  372.         ((eql m1 m2)
  373.          (make-element d1 (lisp:- v1 v2) m1))
  374.         (t (make-element d1 (lisp:- v1 v2) (lisp:gcd m1 m2)))))))
  375.  
  376. (defmethod times ((a GFm-element) (b GFm-element))
  377.   (with-slots ((v1 value) (m1 modulus) (d1 domain)) a
  378.     (with-slots ((v2 value) (m2 modulus) (d2 domain)) b
  379.       (cond ((not (eql d1 d2))
  380.          (error "~S and ~S are not from the same domain" a b))
  381.         ((eql m1 m2)
  382.          (make-element d1 (lisp:* v1 v2) m1))
  383.         (t (make-element d1 (lisp:* v1 v2) (lisp:gcd m1 m2)))))))
  384.  
  385. (defmethod plus ((a GFm-element) (b integer))
  386.   (with-slots ((v1 value) (m1 modulus) (d1 domain)) a
  387.     (make-element d1 (lisp:+ v1 (reduce-modulo-integer b m1)) m1)))
  388.  
  389. (defmethod plus ((a integer) (b GFm-element))
  390.   (with-slots ((v1 value) (m1 modulus) (d1 domain)) b
  391.     (make-element d1 (lisp:+ (reduce-modulo-integer a m1) v1) m1)))
  392.  
  393. (defmethod difference ((a GFm-element) (b integer))
  394.   (with-slots ((v1 value) (m1 modulus) (d1 domain)) a
  395.     (make-element d1 (lisp:- v1 (reduce-modulo-integer b m1)) m1)))
  396.  
  397. (defmethod difference ((a integer) (b GFm-element))
  398.   (with-slots ((v1 value) (m1 modulus) (d1 domain)) b
  399.     (make-element d1 (lisp:- (reduce-modulo-integer a m1) v1) m1)))
  400.  
  401. (defmethod times ((a GFm-element) (b integer))
  402.   (with-slots ((v1 value) (m1 modulus) (d1 domain)) a
  403.     (make-element d1 (lisp:* v1 (reduce-modulo-integer b m1)) m1)))
  404.  
  405. (defmethod times ((a integer) (b GFm-element))
  406.   (with-slots ((v1 value) (m1 modulus) (d1 domain)) b
  407.     (make-element d1 (lisp:* (reduce-modulo-integer a m1) v1) m1)))
  408.  
  409. ;;; Takes the inverse of an integer N mod P.  Solve N*X + P*Y = 1.  N
  410. ;;; is guaranteed to be less than P, since in the case where P is a
  411. ;;; fixnum, N is also assumed to be one.
  412.  
  413. (defmethod recip ((x GFm-element))
  414.   (with-slots (value modulus domain) x
  415.     (make-element domain (reduce-modulo-integer (compute-inverse value modulus)
  416.                         modulus)
  417.           modulus)))
  418.  
  419. (defmethod expt ((x GFm-element) (e integer))
  420.   (with-slots (value modulus domain) x
  421.     (cond ((eql 1 value) x)
  422.       ((lisp:minusp e)
  423.        (error "Raising ~D to a negative power ~D" x e))
  424.       (t (make-element domain (expt-modulo-integer value e modulus)
  425.                modulus)))))
  426.  
  427. (defmethod quotient ((a GFm-element) (b GFm-element)) 
  428.   (with-slots ((v1 value) (m1 modulus) (d1 domain)) a
  429.     (with-slots ((v2 value) (m2 modulus) (d2 domain)) b
  430.       (make-element d1 (lisp:* v1 (compute-inverse v2 m2)) m1))))
  431.   
  432. (defmethod remainder ((a GFm-element) (b GFm-element))
  433.   (error "Computing the remainder of ~D by ~D"
  434.      a b))
  435.  
  436. (defmethod gcd ((a GFm-element) (b GFm-element))
  437.   (with-slots ((v1 value) (m1 modulus) (d1 domain)) a
  438.     (with-slots ((v2 value) (m2 modulus) (d2 domain)) b
  439.       (make-element d1 1 m1))))
  440.  
  441. (defmethod lcm ((a GFm-element) (b GFm-element))
  442.   (with-slots ((v1 value) (m1 modulus) (d1 domain)) a
  443.     (with-slots ((v2 value) (m2 modulus) (d2 domain)) b
  444.       (make-element d1 1 m1))))
  445.  
  446. (defmethod multiplicative-order ((a GFm-element))
  447.   (with-slots (value modulus domain) a
  448.     (cond ((not (eql 1 (lisp:gcd value modulus)))
  449.        *positive-infinity*)
  450.       ((let ((group-order (totient modulus)))
  451.          (do ((factors (factor group-order)
  452.                (rest factors))
  453.           (order group-order))
  454.          ((null factors)
  455.           order)
  456.            (do ((i 0 (lisp:1+ i)))
  457.            ((lisp:= i (cdar factors)))
  458.          (setq order (lisp:/ order (caar factors)))
  459.          (when (not (eql 1 (expt-modulo-integer value order modulus)))
  460.            (setq order (lisp:* order (caar factors)))
  461.            (return t)))))))))
  462.  
  463. (defmethod print-object ((d simple-finite-field) stream)
  464.   #+Genera
  465.   (format stream "~'bGF~(~D)" (characteristic d))
  466.   #-Genera
  467.   (format stream "GF(~D)" (characteristic d)))
  468.  
  469. ;; These are the guys that actually create the finite fields.
  470. (defun make-finite-field* (size)
  471.   (cond ((null size)
  472.      (make-instance 'gfm))
  473.     ((prime? size)
  474.      (make-GFp-domain size 1))
  475.     (t (let* ((s (factor size))
  476.           (char (first (first s)))
  477.           (degree (rest (first s))))
  478.          (if (null (rest s))
  479.          (make-Gfp-domain char degree)
  480.          (error "Finite fields of size ~S=~S don't exist" size s))))))
  481.  
  482.  
  483. (defun make-finite-field (&optional size)
  484.   (add-domain #'false (make-finite-field* size)))
  485.  
  486. ;; This is slightly inefficient, but who cares...  I want to localize
  487. ;; the knowledge of how to create domains in the MAKE-...* functions.
  488. (defun get-finite-field (&optional size)
  489.   (cond ((null size)
  490.      (add-domain (lambda (d) (eql (class-name (class-of d)) 'GFm))
  491.        (make-finite-field* size)))
  492.     ((prime? size)
  493.      (add-domain (lambda (d)
  494.                (and (eql (class-name (class-of d)) 'GFp)
  495.                 (eql (ring-characteristic d) size)))
  496.        (make-finite-field* size)))
  497.      ((null (rest (factor size)))
  498.       (add-domain (lambda (d)
  499.                (and (eql (class-name (class-of d)) 'GF2^n)
  500.                 (eql (lisp:expt (ring-characteristic d)
  501.                             (field-degree d))
  502.                  size)))
  503.        (make-finite-field* size)))
  504.     (t (error "Can't do algebraic extensions yet"))))
  505.